rm(list = ls())
library(DropletUtils)
## Warning: package 'DropletUtils' was built under R version 4.1.2
## Loading required package: SingleCellExperiment
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## Warning: package 'matrixStats' was built under R version 4.1.2
##
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars
## Loading required package: GenomicRanges
## Warning: package 'GenomicRanges' was built under R version 4.1.2
## Loading required package: stats4
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
##
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## Warning: package 'GenomeInfoDb' was built under R version 4.1.2
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians
## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians
library(scater)
## Loading required package: scuttle
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.1.2
library(scran)
## Warning: package 'scran' was built under R version 4.1.2
library(ggplot2)
library(Seurat)
## Warning: package 'Seurat' was built under R version 4.1.2
## Attaching SeuratObject
## Attaching sp
##
## Attaching package: 'Seurat'
## The following object is masked from 'package:SummarizedExperiment':
##
## Assays
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.2
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:Biobase':
##
## combine
## The following objects are masked from 'package:GenomicRanges':
##
## intersect, setdiff, union
## The following object is masked from 'package:GenomeInfoDb':
##
## intersect
## The following objects are masked from 'package:IRanges':
##
## collapse, desc, intersect, setdiff, slice, union
## The following objects are masked from 'package:S4Vectors':
##
## first, intersect, rename, setdiff, setequal, union
## The following objects are masked from 'package:BiocGenerics':
##
## combine, intersect, setdiff, union
## The following object is masked from 'package:matrixStats':
##
## count
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(viridis)
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.1.2
library(EnhancedVolcano)
## Loading required package: ggrepel
## Registered S3 methods overwritten by 'ggalt':
## method from
## grid.draw.absoluteGrob ggplot2
## grobHeight.absoluteGrob ggplot2
## grobWidth.absoluteGrob ggplot2
## grobX.absoluteGrob ggplot2
## grobY.absoluteGrob ggplot2
library(formattable)
##
## Attaching package: 'formattable'
## The following object is masked from 'package:scater':
##
## normalize
## The following object is masked from 'package:BiocGenerics':
##
## normalize
library(pheatmap)
library(ape)
## Warning: package 'ape' was built under R version 4.1.2
library(ComplexHeatmap)
## Loading required package: grid
## ========================================
## ComplexHeatmap version 2.10.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## ! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
## in the original pheatmap() are identically supported in the new function. You
## can still use the original function by explicitly calling pheatmap::pheatmap().
##
## Attaching package: 'ComplexHeatmap'
## The following object is masked from 'package:pheatmap':
##
## pheatmap
library(patchwork)
## Warning: package 'patchwork' was built under R version 4.1.2
##
## Attaching package: 'patchwork'
## The following object is masked from 'package:formattable':
##
## area
library(scales)
## Warning: package 'scales' was built under R version 4.1.2
##
## Attaching package: 'scales'
## The following objects are masked from 'package:formattable':
##
## comma, percent, scientific
## The following object is masked from 'package:viridis':
##
## viridis_pal
library(wesanderson)
library(ggrepel)
library(ggeasy)
library(data.table)
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
## The following object is masked from 'package:SummarizedExperiment':
##
## shift
## The following object is masked from 'package:GenomicRanges':
##
## shift
## The following object is masked from 'package:IRanges':
##
## shift
## The following objects are masked from 'package:S4Vectors':
##
## first, second
library(scuttle)
library(SingleR)
## Warning: package 'SingleR' was built under R version 4.1.2
addTaskCallback(function(...) {set.seed(100);TRUE})
## 2
## 2
##Load Object Object WT and KO values need to be corrected
sce1 <- readRDS("~/BetelLab/WTandReps/Final_Figures/Fayzan.4.SeuratObject.20220420.rds")
cellline <- readRDS("~/BetelLab/WTandReps/Final_Figures/cellline_clean.rds")
fayzan.4 <- sce1
# Take great extra care in counting WT and KO cells!!
# The result should be:
#Always check WT and KO annotations!!
# fayzan.4@meta.data$Sample %>% table
#WT1 WT2 KO1 KO2
#7342 11235 11730 13129
#fayzan.4@meta.data$condition %>% table
#WT KO
#18577 24859
# Before correcting,
fayzan.4@meta.data$sample %>% table
## .
## KO1 KO2 WT1 WT2
## 11235 13129 11730 7342
# KO1 KO2 WT1 WT2
# 11235 13129 11730 7342
fayzan.4@meta.data$condition %>% table
## .
## WT KO
## 19072 24364
# WT KO
# 19072 24364
fayzan.4@meta.data$sampleold = fayzan.4@meta.data$sample
fayzan.4@meta.data$sample = "1"
fayzan.4@meta.data$sample[fayzan.4@meta.data$sampleold == "WT2"] = "WT1"
fayzan.4@meta.data$sample[fayzan.4@meta.data$sampleold == "KO1"] = "WT2"
fayzan.4@meta.data$sample[fayzan.4@meta.data$sampleold == "WT1"] = "KO1"
fayzan.4@meta.data$sample[fayzan.4@meta.data$sampleold == "KO2"] = "KO2"
fayzan.4@meta.data$sample = factor(fayzan.4@meta.data$sample, levels=c("WT1", "WT2", "KO1", "KO2"))
fayzan.4@meta.data$sample %>% table
## .
## WT1 WT2 KO1 KO2
## 7342 11235 11730 13129
# WT1 WT2 KO1 KO2
# 7342 11235 11730 13129
fayzan.4@meta.data$condition = factor(gsub("[12]", "", fayzan.4@meta.data$sample), levels=c("WT", "KO"))
fayzan.4@meta.data$condition %>% table
## .
## WT KO
## 18577 24859
# WT KO
# 18577 24859
fayzan.igraph.csv = fread("Fayzan.4.igraph.cutn.7.20220423.csv")
Idents(fayzan.4) = factor(as.character(fayzan.igraph.csv$igraph_cluster), levels=c("1", "2", "3", "4", "5", "6", "7"))
fayzan.4@meta.data$igraph_cluster = factor(as.character(fayzan.igraph.csv$igraph_cluster), levels=c("1", "2", "3", "4", "5", "6", "7"))
combined_sce1 <- as.SingleCellExperiment(fayzan.4)
combined_sce1$clusters <- fayzan.4@meta.data$igraph_cluster
#saveRDS(combined_sce1, "Final.RDS") #Use condition clusters
#Analysis in Seurat
#sampnames <- combined_sce1$condition
#sampnames2 <- substr(sampnames, 0, 2)
combined_sce2 <- combined_sce1
#gene <- rowData(combined_sce2)$Symbol
#temp <- make.unique(gene, sep ="_")
#rownames(combined_sce2) <- temp
sce2 <- as.Seurat(combined_sce2)
#tmpIdent <- sampnames2
Idents(sce2) <- sce2$clusters
#
sce1 <- sce2
##Figure 4A and Figure 4B
#UMAP
color_samp <- wes_palette("Darjeeling1",3)
wt_col <- color_samp[2]
ko_col <- color_samp[3]
#4A
my_color_palette <- hue_pal()(length(levels(sce1$clusters)))
DimPlot(sce1)
#4B
DimPlot(sce1, group.by = "condition", cols = c(wt_col, ko_col), label = FALSE) + ggtitle("") + theme(legend.text=element_text(size=13))
##Figure 4D
#MAKE HEAT MAP, CLUSTERPLOTS, VLN
markers <- FindAllMarkers(sce1, min.pct = 0.15, min.diff.pct = 0.04, logfc.threshold = .25, test.use = "MAST") #Try with MAST
## Calculating cluster 1
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 2
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 3
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 4
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 5
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 6
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster 7
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
#top20 <- markers %>% dplyr::group_by(markers$cluster) %>% top_n(5, avg_log2FC)
#temp3 <- as.matrix(ScaleData(object = mat, block.size=100))
#png("heatmaptopgroup1month_cluster.png", width = 800, height = 800)
top20 <- markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
#sce1$
sce1 <- ScaleData(object = sce1, block.size=100)
## Centering and scaling data matrix
#png("heatmaptopgroup.png", width = 800, height = 800)
temp3 <- GetAssayData(sce1, slot = "scale.data")
mat<- temp3[unique(c("JUN", top20$gene)), ] %>% as.matrix()
#mat<- t(scale(t(mat)))
heatmapanno <- HeatmapAnnotation(Group = sce1$condition, Cluster = Idents(sce1), col = list(Group=c("WT" = wt_col, "KO" = ko_col), Cluster =c("1" = "#F8766D", "2" = "#C49A00", "3" = "#53B400", "4" = "#00C094", "5" = "#00B6EB", "6" = "#A58AFF", "7" ="#FB61D7")))
col_fun = circlize::colorRamp2(c(-1, 0, 2), c("#140b34" , "#84206b", "#f6d746"))
#Figure 4D
#4D
Heatmap(mat, name = "Scaled Expression",
column_split = Idents(sce1),
#column_title = "Cluster",
cluster_columns = FALSE,
show_column_dend = FALSE,
cluster_column_slices = FALSE,
column_title_gp = gpar(fontsize = 7),
column_gap = unit(.6, "mm"),
cluster_rows = TRUE,
show_row_dend = TRUE,
row_dend_width = unit(20, "mm"),
col = col_fun,
row_names_gp = gpar(fontsize = 8),
column_title_rot = 0,
top_annotation = heatmapanno,
show_column_names = FALSE,
use_raster = TRUE,
raster_quality = 4)
###Figure 4F, exclude cluster 7
sce3 <- subset(sce1, ident = 7, invert = TRUE)
#4F
v1 <- VlnPlot(sce3, "BAX", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##
## This message will be shown once per session.
v2 <- VlnPlot(sce3, "TNFRSF12A", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v3 <- VlnPlot(sce3, "JUN", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v1+v2 + v3
#4F_Alt
v1 <- VlnPlot(sce3, "BAX", split.by = "condition", cols = c(wt_col, ko_col), pt.size = 0)
v2 <- VlnPlot(sce3, "TNFRSF12A", split.by = "condition", cols = c(wt_col, ko_col), pt.size = 0)
v3 <- VlnPlot(sce3, "JUN", split.by = "condition", cols = c(wt_col, ko_col), pt.size = 0)
v1+v2+v3
##Figure 5SG and Figure 4E
Idents(sce1) <-sce1$condition
markers1 <- FindAllMarkers(sce1, min.pct = 0.05, min.diff.pct = 0.02, logfc.threshold = .05, test.use = "MAST")
## Calculating cluster WT
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster KO
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
top <- markers1 %>% dplyr::filter(cluster=="WT")
keyvals <- ifelse(
top$avg_log2FC < -.10, '#F2AD00',
ifelse(top$avg_log2FC > .10, '#00A08A',
'grey50'))
keyvals[is.na(keyvals)] <- 'grey50'
names(keyvals)[keyvals == '#00A08A'] <- 'Wildtype'
names(keyvals)[keyvals == 'grey50'] <- 'NS'
names(keyvals)[keyvals == '#F2AD00'] <- 'Knockout'
apotosis1 <-c("FAS",
"TNFRSF10A",
"TNFRSF10B",
"TNFRSF10C",
"TNFRSF10D",
"TNFRSF11B",
"TNFSF10",
"TNFRSF1A",
"TNFRSF12A",
"FADD",
"CFLAR",
"CASP1",
"CASP2",
"CASP3",
"CASP4",
"CASP6",
"CASP7",
"CASP8",
"CASP9",
"CASP10",
"CASP14",
"NAIP",
"BIRC2",
"BIRC3",
"XIAP",
"BIRC5",
"BIRC6",
"BIRC7",
"BCL2",
"MCL1",
"BCL2L1",
"BCL2L2",
"BCL2A1",
"BCL2L10",
"BAX",
"BAK1",
"BOK",
"BID",
"BCL2L11",
"BMF",
"BAD",
"BIK",
"HRK",
"PMAIP1",
"BNIP3",
"BNIP3L",
"BCL2L14",
"BBC3",
"BCL2L12",
"BCL2L13",
"APAF1",
"CYCS",
"DIABLO",
"HTRA2",
"AIFM1",
"ENDOG",
"CARD8",
"CARD6",
"NOX5",
"TP53",
"CDKN1A",
"CDKN1B",
"CDKN2A",
"CDKN2B")
ofinterest <- c("MKI67",
"TOP2A",
"NES",
"MAP2",
"TH",
"EN1",
"NR4A2",
"HIF1A",
"TP53",
"RELA",
"NFKB1",
"NFKB2",
"JUN",
"FOS", "LMO3", "HES5")
geneofinterest <- c(ofinterest, apotosis1)
#Label genes of interest, change Y axis, Save Size
samplabels1 <- c("PEG10", "DYNC1I2", "JUN",
"NMT1",
"PPDPF",
"PHPT1",
"NDUFA13",
"NDUFC2",
"NDUFA11",
"NME2",
"NDUFA12",
"RPS27L",
"INSM1",
"FABP5",
"HSPA1A",
"WLS",
"MICOS13",
"ASCL1",
"SOX2")
EnhancedVolcano(top,
lab = top$gene,
x = 'avg_log2FC',
y = 'p_val_adj',
title = 'All WT vs KO Top DE Genes',
pCutoff = .0001,
FCcutoff = .10,
pointSize = 2.0,
labSize = 3.3,
col=c('blue', 'green', 'black', 'red3', "yellow", "brown", "grey"),
colAlpha = .8,
selectLab = samplabels1,
colCustom = keyvals,
drawConnectors = TRUE,
widthConnectors = 0.1,
lengthConnectors = unit(0.01, "npc"),
cutoffLineType = "dotted",
cutoffLineCol = "#00000066") + xlim(-.4, .4)
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
sce4 <- subset(sce1, ident == 3 | ident == 5 | ident == 6)
Idents(sce4) <-sce4$condition
markers2 <- FindAllMarkers(sce4, min.pct = 0.01, min.diff.pct = 0.01, logfc.threshold = .05, test.use = "MAST")
## Calculating cluster WT
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster KO
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
top2 <- markers2 %>% dplyr::filter(cluster=="WT")
keyvals <- ifelse(
top2$avg_log2FC > .10 & top2$p_val_adj < .01, '#00A08A',
ifelse(top2$avg_log2FC < -.10 & top2$p_val_adj < .01, '#F2AD00',
'grey50'))
keyvals[is.na(keyvals)] <- 'grey50'
names(keyvals)[keyvals == '#00A08A'] <- 'Wildtype'
names(keyvals)[keyvals == 'grey50'] <- 'NS'
names(keyvals)[keyvals == '#F2AD00'] <- 'Knockout'
samplabels<- c("RPS27L", "PHLDA3", "TUBB2A","TUBB2B","RPL27A","PHPT1","CDKN1A","BAX",
"NMT1",
"BBC3",
"CDKN2B",
"TNFRSF10B")
EnhancedVolcano(top2,
lab = top2$gene,
x = 'avg_log2FC',
y = 'p_val_adj',
title = 'WT vs KO Top DE Genes',
pCutoff = .01,
FCcutoff = .10,
pointSize = 1.0,
labSize = 3.3,
selectLab = samplabels,
col=c('blue', 'green', 'black', 'red3', "yellow", "brown", "grey"),
colAlpha = .8,
colCustom = keyvals,
drawConnectors = TRUE,
widthConnectors = 0.1,
lengthConnectors = unit(0.01, "npc"),
cutoffLineType = "dotted",
cutoffLineCol = "#00000066",
) + xlim(-.75, .75)
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
##Figure 5SE and Figure 5SD #Check Figure Numbering
table2 <- as.matrix(table(sce1@meta.data$condition, sce1@meta.data$cluster))
table3 <- data.table(t(t(table2) / colSums(table2)))
colnames(table3) <- c("Condition", "Cluster", "Fraction")
ggplot(data=table3, aes(x=Cluster, y=Fraction, fill=Condition)) +
geom_bar(stat="identity") + scale_fill_manual(values = c(ko_col, wt_col))
sce1 <- sce2
Idents(sce1) <- sce1$condition
sce1 <- subset(x = sce1, idents = "KO")
th <- length(WhichCells(sce1, slot = 'counts', expression = TH > 0)) / length(WhichCells(sce1, slot = 'counts'))
mki67 <- length(WhichCells(sce1, slot = 'counts', expression = MKI67 > 0)) / length(WhichCells(sce1, slot = 'counts'))
map2 <- length(WhichCells(sce1, slot = 'counts', expression = MAP2 > 0)) / length(WhichCells(sce1, slot = 'counts'))
all_ko <- c(map2, th, mki67)
sce1 <- sce2
Idents(sce1) <- sce1$condition
sce1 <- subset(x = sce1, idents = "WT")
th <- length(WhichCells(sce1, slot = 'counts', expression = TH > 0)) / length(WhichCells(sce1, slot = 'counts'))
mki67 <- length(WhichCells(sce1, slot = 'counts', expression = MKI67 > 0)) / length(WhichCells(sce1, slot = 'counts'))
map2 <- length(WhichCells(sce1, slot = 'counts', expression = MAP2 > 0)) / length(WhichCells(sce1, slot = 'counts'))
all_wt <- c(map2, th, mki67)
temp <- data.table(
Gene = c("MAP2", "TH", "MKI67"),
Group =c("WT", "WT", "WT", "KO", "KO", "KO"),
Values = c(all_wt, all_ko)
)
ggplot(temp, aes(x = Gene, y = Values, fill = Group)) +
geom_bar(stat = "identity", position = "dodge") + ggtitle("Fraction of Cells Expressing Each Gene") + ggeasy::easy_center_title() +
theme(plot.title = element_text(face="bold")) +theme(text=element_text(size=8)) + scale_fill_manual(values = c(ko_col, wt_col))
##Figure 5SF
sce1 <- sce2
sce1.list = list(sce1[, sce1@meta.data$condition == "WT"], sce1[, sce1@meta.data$condition == "KO"])
#MAP2
f1 <-FeaturePlot(sce1.list[[1]], feature="MAP2", order=T, pt.size=1.0) + labs(title="MAP2 (WT)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f2 <-FeaturePlot(sce1.list[[2]], feature="MAP2", order=T, pt.size=1.0) + labs(title="MAP2 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
#PBX1
f3 <-FeaturePlot(sce1.list[[1]], feature="PBX1", order=T, pt.size=1.0) + labs(title="PBX1 (WT)") + scale_colour_gradientn(colours =terrain.colors(10))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f4 <-FeaturePlot(sce1.list[[2]], feature="PBX1", order=T, pt.size=1.0) + labs(title="PBX1 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
#MKI67
f5 <-FeaturePlot(sce1.list[[1]], feature="MKI67", order=T, pt.size=1.0) + labs(title="MKI67 (WT)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f6 <-FeaturePlot(sce1.list[[2]], feature="MKI67", order=T, pt.size=1.0) + labs(title="MKI67 (KO)") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
(f1/f2) | (f3/f4) | (f5/f6)
##Figure 4C and Figure 4H La Manno Annotations
combined_sce1 <- as.SingleCellExperiment(sce1)
combined_sce1$clusters <- sce1$clusters
combined_sce_lamanno <- scRNAseq::LaMannoBrainData("human-embryo")
## snapshotDate(): 2021-10-19
## see ?scRNAseq and browseVignettes('scRNAseq') for documentation
## loading from cache
## see ?scRNAseq and browseVignettes('scRNAseq') for documentation
## loading from cache
# SingleR() expects reference datasets to be normalized and log-transformed.
combined_sce_lamanno <- logNormCounts(combined_sce_lamanno)
predictions2 <- SingleR(test=combined_sce1, assay.type.test=2,
ref=combined_sce_lamanno, labels =combined_sce_lamanno$Cell_type, de.method="wilcox", clusters=combined_sce1$clusters)
predictions2$pruned.labels #Didn't assign
## [1] "hNbM" "hNbM" "hProgFPL" "hNbM" "hProgFPL" "hProgFPL" "hPeric"
sce5 <- as.Seurat(combined_sce1)
Idents(sce5) <- sce5$clusters
sce5 <- RenameIdents(object = sce5, `1` = "hNbM", `2` = "hNbM", `3` = "hProgFPL", '4'= "hNbM", "5" = "hProgFPL", "6" = "hProgFPL", "7" = "hPeric")
DimPlot(sce5)
#4H by celltype
combined_sce_lamanno1 <- subset(combined_sce_lamanno , ,Cell_type == c('hNbM') | Cell_type == c('hNbML1') | Cell_type == c('hProgFPL'))
predictions3 <- SingleR(test=combined_sce1, assay.type.test=2,
ref=combined_sce_lamanno1, labels =combined_sce_lamanno1$Cell_type, de.method="wilcox")
combined_sce1$pruned.labels <- predictions3$pruned.labels
sce6 <- as.Seurat(combined_sce1)
#cols2 <- wes_palette("Rushmore1",3)
v1 <- VlnPlot(sce6, "BAX", group.by = "pruned.labels", cols = c(wt_col, ko_col), pt.size = 0, split.by = "condition")
v2 <- VlnPlot(sce6, "TNFRSF12A", group.by = "pruned.labels", cols = c(wt_col, ko_col), pt.size = 0, split.by = "condition")
v3 <- VlnPlot(sce6, "JUN", group.by = "pruned.labels", cols = c(wt_col, ko_col), pt.size = 0, split.by = "condition")
v1 + v2 + v3
##Figure 4G, Use Cellline Object
#CellLine
#Load Cellline
#cellline <- readRDS("cellline_clean.rds")
#HES5
f1 <- FeaturePlot(sce1, feature="HES5", order=T, pt.size=1.0) + labs(title="HES5") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f2 <- FeaturePlot(cellline, feature="HES5", order=T, pt.size=1.0) + labs(title="HES5") + scale_colour_gradientn(colours = terrain.colors(10))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
f1|f2
##Figure ALT
#ALTERNATE FIGURE REQUESTS
sce4 <- subset(sce1, ident == 1 |ident == 2 | ident == 4)
Idents(sce4) <-sce4$condition
markers2 <- FindAllMarkers(sce4, min.pct = 0.01, min.diff.pct = 0.01, logfc.threshold = .05, test.use = "MAST")
## Calculating cluster WT
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
## Calculating cluster KO
##
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
##
## Done!
top2 <- markers2 %>% dplyr::filter(cluster=="WT")
keyvals <- ifelse(
top2$avg_log2FC >= .10 & top2$p_val_adj < .01, '#00A08A',
ifelse(top2$avg_log2FC <= -.10 & top2$p_val_adj < .01, '#F2AD00',
'grey50'))
keyvals[is.na(keyvals)] <- 'grey50'
names(keyvals)[keyvals == '#00A08A'] <- 'Wildtype'
names(keyvals)[keyvals == 'grey50'] <- 'NS'
names(keyvals)[keyvals == '#F2AD00'] <- 'Knockout'
samplabels<- c("PEG10",
"NMT1",
"HSPA8", "RPS27L",
"NDUFA12",
"DYNLT1",
"PDCD5",
"PHPT1",
"KRTCAP2", "SYT14",
"SOX2",
"HSPA1B",
"HSPA1A","HSPA2A","DNAJB1","HSPA8")
EnhancedVolcano(top2,
lab = top2$gene,
x = 'avg_log2FC',
y = 'p_val_adj',
title = 'WT vs KO Top DE Genes',
pCutoff = .01,
FCcutoff = .10,
pointSize = 1.0,
labSize = 3.3,
selectLab = samplabels,
col=c('blue', 'green', 'black', 'red3', "yellow", "brown", "grey"),
colAlpha = .8,
colCustom = keyvals,
drawConnectors = TRUE,
widthConnectors = 0.1,
lengthConnectors = unit(0.01, "npc"),
cutoffLineType = "dotted",
cutoffLineCol = "#00000066",
#titleLabSize = 10,
#subtitleLabSize = 8,
#captionLabSize = 8
) + xlim(-.8, .8)
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
##### ALTERNATE VIOLIN PLOTS
sce3 <- subset(sce1, ident = 7, invert = TRUE)
#ATF3, ATF4, JUN, FOS, HSPD1, HSPA1A, HSPA2A, and DNAJB1 HSPA8 DNAJB2 ATF2
v1 <- VlnPlot(sce3, "ATF3", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v2 <- VlnPlot(sce3, "ATF4", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v3 <- VlnPlot(sce3, "JUN", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v4 <- VlnPlot(sce3, "FOS", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v5 <- VlnPlot(sce3, "HSPD1", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v6 <- VlnPlot(sce3, "HSPA1A", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v8 <- VlnPlot(sce3, "DNAJB1", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v9 <- VlnPlot(sce3, "HSPA8", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v10 <- VlnPlot(sce3, "DNAJB2", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v11 <- VlnPlot(sce3, "ATF2", split.by = "condition", split.plot = TRUE, cols = c(wt_col, ko_col), pt.size = 0)
v8+v9